// -*-C++-*-------------------------------------------------------------------
//
// author: Silvio Tosatto
//
// date: 03/03
//
// descriptipton: defines common statistical operations, e.g. AVG, SD
//
// modified by: Luca Stoppa on 2/4/2008 for random numbers generator.
//
// ---------------------------------------------------------------------------

#ifndef __STAT_TOOLS_H__
#define __STAT_TOOLS_H__

#include <math.h>
#include <vector>
#include <Debug.h>
#include <algorithm>

template <class T> 
inline 
double minimum(vector<T> data)
{
  T min = data[0];
  for (unsigned int i = 1; i < data.size(); i++)
    if (min > data[i])
      min = data[i];
  return min;
}

template <class T> 
inline 
double maximum(vector<T> data)
{
  T max = data[0];
  for (unsigned int i = 1; i < data.size(); i++)
    if (max < data[i])
      max = data[i];
  return max;
}

template <class T> 
inline 
double average(vector<T> data)
{
  T avg = 0;
  for (unsigned int i = 0; i < data.size(); i++)
      avg += data[i];
  avg = avg / data.size();
  return avg;
}

template <class T> 
inline 
double standardDeviation(vector<T> data, T avg = -9999)
{
  if (avg == -9999)
    avg = average(data);

  T sd = 0;
  for (unsigned int i = 0; i < data.size(); i++)
      sd += (data[i] - avg) * (data[i] - avg);
  sd = sqrt(sd / data.size());
  return sd;
}

template <class T> 
inline 
double standardError(vector<T> data, T sd = -9999)
{
  if (sd == -9999)
    sd = standardDeviation(data);

  return (sd / sqrt(data.size()));
}


template <class T> 
inline 
vector<double> Zscore(vector<T> data, T avg = -9999, T sd = -9999)
{
  if (avg == -9999)
    avg = average(data);
  if (sd == -9999)
    sd = standardDeviation(data, avg);

  vector<double> tmp;
  tmp.reserve(data.size());

  for (unsigned int i = 0; i < data.size(); i++)
    {
      tmp.push_back( (data[i]-avg) / sd  );
    }

  return tmp;
}


template <class T> 
inline 
vector<double> windowAverage(vector<T> data, unsigned int w)
{
  vector<double> tmp;
  tmp.reserve(data.size());

  for (int i = 0; i < static_cast<int>(data.size()); i++)
    {
      double val = 0.0;

      for (unsigned int j = i - w; j <= i + w; j++)
	{
	  if (j < 0)
	    val = data[0] + val;
	  else if (j >= data.size())
	    val = data[data.size()-1] + val;
	  else
	    val = data[j] + val;
	}

      tmp.push_back( val / (2 * w + 1)  );
    }

  return tmp;
}

template <class T> 
inline 
double pearsonCorrelation(vector<T> data, vector<T> data2)
{
  T pearson = 0;
  long double agv1 = average(data);
  long double agv2 = average(data2);
  long double numerator = 0;
  for ( unsigned int i = 0; i < data.size(); i++)
    {
      numerator = numerator + ((data[i] - agv1 ) * (data2[i] - agv2));
    }
  long double tmpdenom1 = 0;
  long double denom1 = 0;
  for ( unsigned int i = 0; i < data.size(); i++ )
    {
      tmpdenom1 = tmpdenom1 + ( (data[i] - agv1) * (data[i] - agv1) );
    }
  denom1 = sqrt(tmpdenom1);

  long double tmpdenom2 = 0;
  long double denom2 = 0;
  for ( unsigned int i = 0; i < data2.size(); i++ )
    {
      tmpdenom2 = tmpdenom2 + ( (data2[i] - agv2) * (data2[i] - agv2) );
    }
  denom2 = sqrt(tmpdenom2);
  
  pearson = ( numerator / ( denom1 * denom2 ) );
  return pearson;
}

template <class T>
inline
double sRankHelper ( T n, vector<T> data )
{
  double result;
  for (unsigned int i = 0; i < data.size(); i++ )
    {
      if ( n == data[i] )
	{
	  result = static_cast<double>(i);
	  return result;
	}
    }
  return -10000.00;
}


template <class T>
inline 
double spearmanCorrelation(vector<T> data, vector<T> data2)
{
  vector<T> ordData; 
  vector<T> ordData2;
  
  for ( unsigned int i = 0; i < data.size(); i++ ) 
    {
      ordData.push_back(data[i]);
      ordData2.push_back(data2[i]);
    }
  
  sort(&ordData[0], &ordData[ordData.size()-1]);
  sort(&ordData2[0], &ordData2[ordData2.size()-1]);
  
  T spearman = 0;
  long double D = 0;
  for ( unsigned int i = 0; i < data.size(); i++ )
    D = ( ( sRankHelper(data[i], ordData) - 
	    ( sRankHelper(data2[i], ordData2) )) * 
	  ( sRankHelper(data[i], ordData) - 
	    ( sRankHelper(data2[i], ordData2) ))); 
  
  spearman = 1 - ( ( 6*D ) / ( ( data.size()*data.size()*data.size() ) 
			       - data.size() ) );
  return spearman;
}


// ---------------------------------------------------------------------------
// NAME: getRandomNumber()
// DESC: given an interval [a,b], we generate random number inside it.
// NOTE: before calling this function, please initialize a random seed by
//       using srand( (unsigned)time(0) ) ONCE in your program.
// AUTH: Luca Stoppa
// DATE: 3/4/2008
// ---------------------------------------------------------------------------
double getRandomNumber( double low, double high );

// ---------------------------------------------------------------------------
// NAME: getGaussianRandomNumber()
// DESC: given an interval [a,b], we generate a random number inside it based
//       on a gaussian distribution of probability.
// NOTE: before calling this function, please initialize a random seed by
//       using srand( (unsigned)time(0) ) ONCE in your program.
//       Be aware on using as seed the value 12: this is an experimental
//       result that gets almost the best distribution of probability.
// AUTH: Luca Stoppa
// DATE: 3/4/2008
// ---------------------------------------------------------------------------
double getGaussianRandomNumber( double low, double high, double seed=12 );


#endif // __STAT_TOOLS_H__
